GAMMs analyses Trombone - Tongan vs English (NZE) -> Hypothesis 1, Prediction 1b

This notebook provides the analysis regarding Hypothesis 1, Prediction 1b for the article: Heyne, M., Derrick, D., and Al-Tamimi, J. (under review). “Native language influence on brass instrument performance: An application of generalized additive mixed models (GAMMs) to midsagittal ultrasound images of the tongue”. Frontiers Research Topic: Models and Theories of Speech Production. Ed. Adamantios Gafos & Pascal van Lieshout.

# specify directory to save models and summaries
output_dir = "updated_models"
# specify whether to run models -> if set to false script will attempt to load saved models from output_dir
run_models = FALSE

1 Loading packages

load_packages = c("parallel","mgcv", "itsadug", "rlist", "ggplot2", "plotly", "dplyr")
for(pkg in load_packages){
  eval(bquote(library(.(pkg))))
  if (paste0("package:", pkg) %in% search()){
    cat(paste0("Successfully loaded the ", pkg, " package.\n"))
  }else{
    install.packages(pkg)
    eval(bquote(library(.(pkg))))
    if (paste0("package:", pkg) %in% search()){
      cat(paste0("Successfully loaded the ", pkg, " package.\n"))
    }
  }
}

Successfully loaded the parallel package.
Successfully loaded the mgcv package.
Successfully loaded the itsadug package.
Successfully loaded the rlist package.
Successfully loaded the ggplot2 package.
Successfully loaded the plotly package.
Successfully loaded the dplyr package.

rm(load_packages, pkg)
# detect number of cores available for model calculations
ncores = detectCores()
cat(paste0("Number of cores available for model calculations set to ", ncores, "."))

Number of cores available for model calculations set to 8.

1.0.1 plotly_model_outputs function (Matthias Heyne, 2019)

# plot multiple GAM model outputs
plotly_model_outputs <- function(model, changing_cond, changing_var, constant_cond1, constant_var1, values, constant_cond2=NULL, constant_var2=NULL, print=TRUE){
  if(length(constant_var1)>1 | length(constant_var2)>1){
    print("Error: Constant variables can only have length 1.")
  }else{
    if(!is.null(constant_cond2) && !is.null(constant_var2) && length(changing_var)==2){
      # works for models Notes.gam...
      cond_p1 = capture.output(cat(paste0("list(", changing_cond, "='", changing_var[1], "', ", constant_cond1, "='", constant_var1, "', ", constant_cond2, "='", constant_var2, "')")))
      p1=plot_smooth(x=get(model), view=values, cond = eval(parse(text=cond_p1)))
      cond_p2 = capture.output(cat(paste0("list(", changing_cond, "='", changing_var[2], "', ", constant_cond1, "='", constant_var1, "', ", constant_cond2, "='", constant_var2, "')")))
      p2=plot_smooth(x=get(model), view=values, cond = eval(parse(text=cond_p2)))
      
      # set Rho max to the max of the fit + half the difference between max of the fit and the upper limit
      max_ul = max(p1$fv$ul, p2$fv$ul)
      max_fit = max(p1$fv$fit, p2$fv$fit)
      maximum=max_fit+((max_ul-max_fit)/2)
      
      # plot in polar coordinates
      p=plot_ly(type='scatterpolar', mode='lines') %>% 
        add_trace(theta=p1$fv$theta_uncut_z*180/pi, r=p1$fv$fit, line=list(color="blue", width=2.5), name=changing_var[1]) %>% 
        add_trace(theta=p1$fv$theta_uncut_z*180/pi, r=p1$fv$ul, line=list(color="blue", dash="dot", width=0.5), showlegend=FALSE) %>% 
        add_trace(theta=p1$fv$theta_uncut_z*180/pi, r=p1$fv$ll, line=list(color="blue", dash="dot", width=0.5), showlegend=FALSE) %>% 
        add_trace(theta=p2$fv$theta_uncut_z*180/pi, r=p2$fv$fit, line=list(color="red", width=2.5), name=changing_var[2]) %>% 
        add_trace(theta=p2$fv$theta_uncut_z*180/pi, r=p2$fv$ul, line=list(color="red", dash="dot", width=0.5), showlegend=FALSE) %>% 
        add_trace(theta=p2$fv$theta_uncut_z*180/pi, r=p2$fv$ll, line=list(color="red", dash="dot", width=0.5), showlegend=FALSE) %>%
        layout(polar=list(sector=c(20,160), radialaxis=list(angle=90, range=c(0,maximum)),
                          angularaxis=list(thetaunit='radians', direction="clockwise", rotation=0)), 
               title=paste0("GAM smooths @", constant_var1, " & ", constant_var2))
      p
    }else if(is.null(constant_cond2) && is.null(constant_var2) && length(changing_var)==2){
      # no specific case yet
      cond_p1 = capture.output(cat(paste0("list(", changing_cond, "='", changing_var[1], "', ", constant_cond1, "='", constant_var1, "')")))
      p1=plot_smooth(x=get(model), view=values, cond = eval(parse(text=cond_p1)))
      cond_p2 = capture.output(cat(paste0("list(", changing_cond, "='", changing_var[2], "', ", constant_cond1, "='", constant_var1, "')")))
      p2=plot_smooth(x=get(model), view=values, cond = eval(parse(text=cond_p2)))
      
      # set Rho max to the max of the fit + half the difference between max of the fit and the upper limit
      max_ul = max(p1$fv$ul, p2$fv$ul)
      max_fit = max(p1$fv$fit, p2$fv$fit)
      maximum=max_fit+((max_ul-max_fit)/2)
      
      # plot in polar coordinates
      p=plot_ly(type='scatterpolar', mode='lines') %>% 
        add_trace(theta=p1$fv$theta_uncut_z*180/pi, r=p1$fv$fit, line=list(color="blue", width=2.5), name=changing_var[1]) %>% 
        add_trace(theta=p1$fv$theta_uncut_z*180/pi, r=p1$fv$ul, line=list(color="blue", dash="dot", width=0.5), showlegend=FALSE) %>% 
        add_trace(theta=p1$fv$theta_uncut_z*180/pi, r=p1$fv$ll, line=list(color="blue", dash="dot", width=0.5), showlegend=FALSE) %>% 
        add_trace(theta=p2$fv$theta_uncut_z*180/pi, r=p2$fv$fit, line=list(color="red", width=2.5), name=changing_var[2]) %>% 
        add_trace(theta=p2$fv$theta_uncut_z*180/pi, r=p2$fv$ul, line=list(color="red", dash="dot", width=0.5), showlegend=FALSE) %>% 
        add_trace(theta=p2$fv$theta_uncut_z*180/pi, r=p2$fv$ll, line=list(color="red", dash="dot", width=0.5), showlegend=FALSE) %>%
        layout(polar=list(sector=c(20,160), radialaxis=list(angle=90, range=c(0,maximum)),
                          angularaxis=list(thetaunit='radians', direction="clockwise", rotation=0)), 
               title=paste0("GAM smooths @", constant_var1, " & ", constant_var2))
      p
    }else if(is.null(constant_cond2) && is.null(constant_var2) && length(changing_var)>2){
      # works for models NZE.gam... or Tongan.gam...
      for (i in 1:length(changing_var)){
        if(changing_cond == constant_cond1){
          cond_p1 = capture.output(cat(paste0("list(", constant_cond1, "='", constant_var1, "')")))
          p1=plot_smooth(x=get(model), view=values, cond = eval(parse(text=cond_p1)))
          # exception for KIT (='\\\\') when not using IPA symbols
          if (changing_var[i]!="\\\\"){
            cond_p2 = capture.output(cat(paste0("list(", changing_cond, "='", changing_var[i], "')")))
            p2=plot_smooth(x=get(model), view=values, cond = eval(parse(text=cond_p2)))
          }else{
            p2=plot_smooth(x=get(model), view=values, cond = list(token.ord='\\\\'))
          }
          
          # set Rho max to the max of the fit + half the difference between max of the fit and the upper limit
          max_ul = max(p1$fv$ul, p2$fv$ul)
          max_fit = max(p1$fv$fit, p2$fv$fit)
          maximum=max_fit+((max_ul-max_fit)/2)
          # plot in polar coordinates
          p=plot_ly(type='scatterpolar', mode='lines') %>% 
            add_trace(theta=p1$fv$theta_uncut_z*180/pi, r=p1$fv$fit, line=list(color="blue", width=2.5), name=constant_var1) %>% 
            add_trace(theta=p1$fv$theta_uncut_z*180/pi, r=p1$fv$ul, line=list(color="blue", dash="dot", width=0.5), showlegend=FALSE) %>% 
            add_trace(theta=p1$fv$theta_uncut_z*180/pi, r=p1$fv$ll, line=list(color="blue", dash="dot", width=0.5), showlegend=FALSE) %>% 
            add_trace(theta=p2$fv$theta_uncut_z*180/pi, r=p2$fv$fit, line=list(color="red", width=2.5), name=changing_var[i]) %>% 
            add_trace(theta=p2$fv$theta_uncut_z*180/pi, r=p2$fv$ul, line=list(color="red", dash="dot", width=0.5), showlegend=FALSE) %>% 
            add_trace(theta=p2$fv$theta_uncut_z*180/pi, r=p2$fv$ll, line=list(color="red", dash="dot", width=0.5), showlegend=FALSE) %>%
            layout(polar=list(sector=c(20,160), radialaxis=list(angle=90, range=c(0,maximum)),
                              angularaxis=list(thetaunit='radians', direction="clockwise", rotation=0)), 
                   title=paste0("GAM smooths ",constant_var1," vs ", changing_var[i]))
          Sys.sleep(0)
          print(p)
        }
      }
    }
  }
}

1.0.2 plotly_sig_diff_notes function (Matthias Heyne, 2019) -> with updated legend position

# plot smooths with shading to indicate significant differences (Matthias Heyne, 2019)
plotly_smooths_w_sig_diff <- function(model, condition, var1, var2, values, language, print=TRUE){
  # get intervals of significant differences by running plot_diff
  # unfortunately setting plot=FALSE doesn't work as intervals of significant difference are not displayed!
  # hardcoded condition
  # output = capture.output(plot_diff(get(model), view=values, 
  #                                   comp=list(tokenPooled.ord=c(var1, var2))))
  # output = capture.output(plot_diff(get(model), view=values, comp=list(langNoteInt.ord=c(paste0("Tongan.", note, ".", intensity), paste0("NZE.", note, ".", intensity)))))
  names_smooths=list()
  if (condition=="tokenPooled.ord" && length(language)==1){
    output_comp = capture.output(cat(paste0("list(", condition, "=c(var1, var2))")))
    names_smooths[1]=var1
    names_smooths[2]=var2
    plot_title = paste0("GAM smooths ", language, " ", var1, " vs ", var2)
  }else if (condition=="langNoteInt.ord" && length(language)==2){
    output_comp = capture.output(cat(paste0("list(", condition, "=c('", language[1], ".", var1, ".", var2, 
                                            "', '", language[2], ".", var1, ".", var2, "'))")))
    names_smooths[1]=paste0(language[1], ".", var1, ".", var2)
    names_smooths[2]=paste0(language[2], ".", var1, ".", var2)
    plot_title = paste0("GAM smooths ", language[1], ".", var1, ".", var2, " vs ", language[2], ".", var1, ".", var2)
  }else if (condition=="native_lg.ord"){
    output_comp = capture.output(cat(paste0("list(", condition, "=c(var1, var2))")))
    names_smooths[1]=var1
    names_smooths[2]=var2
    plot_title = paste0("GAM smooths ", language, " ", var1, " vs ", var2)
  }
  
  # output_comp = capture.output(cat(paste0("list(", condition, "=c(var1, var2))")))
  output = capture.output(plot_diff(get(model), view=values, comp=eval(parse(text=output_comp))))
  # no significant difference
  if ((length(language)==1 && length(output)==7) | (length(language)==2 && length(output)==6)){
    cat(paste0("Smooths for ", var1, " & ", var2, " are not significantly different.\n"))
    dat1 = NA
    assign(paste0("int_sig_diff_", var1, "_", var2), dat1, envir = .GlobalEnv)
    rm(dat1)
    # run plot_smooth to grab data for polar plots
    p1 = plot_smooth(x=get(model), view=values, cond=list(tokenPooled.ord=var1, tokenPooled.ord=var2))
    p2 = plot_smooth(x=get(model), view=values, cond=list(tokenPooled.ord=var2, tokenPooled.ord=var1))
    # set Rho max to the max of the fit + half the difference between max of the fit and the upper limit
    maximum=max(p1$fv$fit, p2$fv$fit)+((max(p1$fv$ul, p2$fv$ul)-max(p1$fv$fit, p2$fv$fit))/2)
    # plot in polar coordinates
    p=plot_ly(type='scatterpolar', mode='lines') %>% 
      add_trace(theta=p1$fv$theta_uncut_z*180/pi, r=p1$fv$fit, line=list(color="blue", width=2.5), name=names_smooths[1]) %>% 
      add_trace(theta=p1$fv$theta_uncut_z*180/pi, r=p1$fv$ul, line=list(color="blue", dash="dot", width=0.5), showlegend=FALSE) %>% 
      add_trace(theta=p1$fv$theta_uncut_z*180/pi, r=p1$fv$ll, line=list(color="blue", dash="dot", width=0.5), showlegend=FALSE) %>% 
      add_trace(theta=p2$fv$theta_uncut_z*180/pi, r=p2$fv$fit, line=list(color="red", width=2.5), name=names_smooths[2]) %>% 
      add_trace(theta=p2$fv$theta_uncut_z*180/pi, r=p2$fv$ul, line=list(color="red", dash="dot", width=0.5), showlegend=FALSE) %>% 
      add_trace(theta=p2$fv$theta_uncut_z*180/pi, r=p2$fv$ll, line=list(color="red", dash="dot", width=0.5), showlegend=FALSE) %>%
      layout(polar=list(sector=c(20,160), radialaxis=list(angle=90, range=c(0,maximum)),
                        angularaxis=list(thetaunit='radians', direction="clockwise", rotation=0)), 
             title=plot_title, legend=list(orientation="h", xanchor="center", x=0.5))
    Sys.sleep(0)
    print(p)
  # there are differences...
  }else{
    # grab intervals of significant differences from output
    if (length(language)==1 && length(output)>=8){
      sig_diff1 = c(as.double(unlist(strsplit(unlist(strsplit(output[8], " "))[1], "\t"))[2]), as.double(unlist(strsplit(output[8], " "))[3]))
    }else if (length(language)==2 && length(output)>=7){
      sig_diff1 = c(as.double(unlist(strsplit(unlist(strsplit(output[7], " "))[1], "\t"))[2]), as.double(unlist(strsplit(output[7], " "))[3]))
    }
    if (length(language)==1 && length(output)>=9){
      sig_diff2 = c(as.double(unlist(strsplit(unlist(strsplit(output[9], " "))[1], "\t"))[2]), as.double(unlist(strsplit(output[9], " "))[3]))
    }else if (length(language)==2 && length(output)>=8){
      sig_diff2 = c(as.double(unlist(strsplit(unlist(strsplit(output[8], " "))[1], "\t"))[2]), as.double(unlist(strsplit(output[8], " "))[3]))
    }
    if (length(language)==1 && length(output)>=10){
      sig_diff3 = c(as.double(unlist(strsplit(unlist(strsplit(output[10], " "))[1], "\t"))[2]), as.double(unlist(strsplit(output[10], " "))[3]))
    }else if (length(language)==2 && length(output)>=9){
      sig_diff3 = c(as.double(unlist(strsplit(unlist(strsplit(output[9], " "))[1], "\t"))[2]), as.double(unlist(strsplit(output[9], " "))[3]))
    }
    if (length(language)==1 && length(output)>=11){
      sig_diff4 = c(as.double(unlist(strsplit(unlist(strsplit(output[11], " "))[1], "\t"))[2]), as.double(unlist(strsplit(output[11], " "))[3]))
    }else if (length(language)==2 && length(output)>=10){
      sig_diff4 = c(as.double(unlist(strsplit(unlist(strsplit(output[10], " "))[1], "\t"))[2]), as.double(unlist(strsplit(output[10], " "))[3]))
    }
    # write intervals of significant difference to variable
    if ((length(language)==1 && length(output)>=11) | (length(language)==2 && length(output)>=10)){
      dat1 = c(sig_diff1, sig_diff2, sig_diff3, sig_diff4)
    }else if ((length(language)==1 && length(output)>=10) | (length(language)==2 && length(output)>=9)){
      dat1 = c(sig_diff1, sig_diff2, sig_diff3)
    }else if ((length(language)==1 && length(output)>=9) | (length(language)==2 && length(output)>=8)){
      dat1 = c(sig_diff1, sig_diff2)
    }else{
      dat1 = sig_diff1
    }
    assign(paste0("int_sig_diff_", var1, "_", var2), dat1, envir = .GlobalEnv)
    rm(dat1, output)
    # run plot_smooth to grab data for polar plots
    if (condition=="tokenPooled.ord" && length(language)==1){
      cond_p1 = capture.output(cat(paste0("list(", condition, "=var1)")))
      cond_p2 = capture.output(cat(paste0("list(", condition, "=var2)")))
    }else if (condition=="langNoteInt.ord" && length(language)==2){
      cond_p1 = capture.output(cat(paste0("list(", condition, "='", language[1], ".", var1, ".", var2, 
                                "', ", condition, "='", language[2], ".", var1, ".", var2, "')")))
      cond_p2 = capture.output(cat(paste0("list(", condition, "='", language[2], ".", var1, ".", var2, 
                                          "', ", condition, "='", language[1], ".", var1, ".", var2, "')")))
    }else if (condition=="native_lg.ord"){
      cond_p1 = capture.output(cat(paste0("list(", condition, "=var1)")))
      cond_p2 = capture.output(cat(paste0("list(", condition, "=var2)")))
    }
    # hardcoded condition
    p1 = plot_smooth(x=get(model), view=values, cond=eval(parse(text=cond_p1)))
    p2 = plot_smooth(x=get(model), view=values, cond=eval(parse(text=cond_p2)))
    # set Rho max to the max of the fit + half the difference between max of the fit and the upper limit
    maximum=max(p1$fv$fit, p2$fv$fit)+((max(p1$fv$ul, p2$fv$ul)-max(p1$fv$fit, p2$fv$fit))/2)
    # plot in polar coordinates
    if (exists("sig_diff4")){
      p=plot_ly(type='scatterpolar', mode='lines') %>% 
        add_trace(theta=p1$fv$theta_uncut_z*180/pi, r=p1$fv$fit, line=list(color="blue", width=2.5), name=names_smooths[1]) %>% 
        add_trace(theta=p1$fv$theta_uncut_z*180/pi, r=p1$fv$ul, line=list(color="blue", dash="dot", width=0.5), showlegend=FALSE) %>% 
        add_trace(theta=p1$fv$theta_uncut_z*180/pi, r=p1$fv$ll, line=list(color="blue", dash="dot", width=0.5), showlegend=FALSE) %>% 
        add_trace(theta=p2$fv$theta_uncut_z*180/pi, r=p2$fv$fit, line=list(color="red", width=2.5), name=names_smooths[2]) %>% 
        add_trace(theta=p2$fv$theta_uncut_z*180/pi, r=p2$fv$ul, line=list(color="red", dash="dot", width=0.5), showlegend=FALSE) %>% 
        add_trace(theta=p2$fv$theta_uncut_z*180/pi, r=p2$fv$ll, line=list(color="red", dash="dot", width=0.5), showlegend=FALSE) %>%
        add_trace(theta=seq(sig_diff1[1]*180/pi, sig_diff1[2]*180/pi, length.out=20), 
                  r=c(0, rep(maximum, 18), 0), line=list(color="black", width=0.5), fill="toself", fillcolor=rgb(0,0,0,max=255,alpha=25), showlegend=FALSE) %>%
        add_trace(theta=seq(sig_diff2[1]*180/pi, sig_diff2[2]*180/pi, length.out=20), 
                  r=c(0, rep(maximum, 18), 0), line=list(color="black", width=0.5), fill="toself", fillcolor=rgb(0,0,0,max=255,alpha=25), showlegend=FALSE) %>%
        add_trace(theta=seq(sig_diff3[1]*180/pi, sig_diff3[2]*180/pi, length.out=20), 
                  r=c(0, rep(maximum, 18), 0), line=list(color="black", width=0.5), fill="toself", fillcolor=rgb(0,0,0,max=255,alpha=25), showlegend=FALSE) %>%
        add_trace(theta=seq(sig_diff4[1]*180/pi, sig_diff4[2]*180/pi, length.out=20), 
                  r=c(0, rep(maximum, 18), 0), line=list(color="black", width=0.5), fill="toself", fillcolor=rgb(0,0,0,max=255,alpha=25), showlegend=FALSE) %>%
        layout(polar=list(sector=c(20,160), radialaxis=list(angle=90, range=c(0,maximum)),
                          angularaxis=list(thetaunit='radians', direction="clockwise", rotation=0)), 
               title=plot_title, legend=list(orientation="h", xanchor="center", x=0.5))
      Sys.sleep(0)
      print(p)
    }else if (exists("sig_diff3")){
      p=plot_ly(type='scatterpolar', mode='lines') %>% 
        add_trace(theta=p1$fv$theta_uncut_z*180/pi, r=p1$fv$fit, line=list(color="blue", width=2.5), name=names_smooths[1]) %>% 
        add_trace(theta=p1$fv$theta_uncut_z*180/pi, r=p1$fv$ul, line=list(color="blue", dash="dot", width=0.5), showlegend=FALSE) %>% 
        add_trace(theta=p1$fv$theta_uncut_z*180/pi, r=p1$fv$ll, line=list(color="blue", dash="dot", width=0.5), showlegend=FALSE) %>% 
        add_trace(theta=p2$fv$theta_uncut_z*180/pi, r=p2$fv$fit, line=list(color="red", width=2.5), name=names_smooths[2]) %>% 
        add_trace(theta=p2$fv$theta_uncut_z*180/pi, r=p2$fv$ul, line=list(color="red", dash="dot", width=0.5), showlegend=FALSE) %>% 
        add_trace(theta=p2$fv$theta_uncut_z*180/pi, r=p2$fv$ll, line=list(color="red", dash="dot", width=0.5), showlegend=FALSE) %>%
        add_trace(theta=seq(sig_diff1[1]*180/pi, sig_diff1[2]*180/pi, length.out=20), 
                  r=c(0, rep(maximum, 18), 0), line=list(color="black", width=0.5), fill="toself", fillcolor=rgb(0,0,0,max=255,alpha=25), showlegend=FALSE) %>%
        add_trace(theta=seq(sig_diff2[1]*180/pi, sig_diff2[2]*180/pi, length.out=20), 
                  r=c(0, rep(maximum, 18), 0), line=list(color="black", width=0.5), fill="toself", fillcolor=rgb(0,0,0,max=255,alpha=25), showlegend=FALSE) %>%
        add_trace(theta=seq(sig_diff3[1]*180/pi, sig_diff3[2]*180/pi, length.out=20), 
                  r=c(0, rep(maximum, 18), 0), line=list(color="black", width=0.5), fill="toself", fillcolor=rgb(0,0,0,max=255,alpha=25), showlegend=FALSE) %>%
        layout(polar=list(sector=c(20,160), radialaxis=list(angle=90, range=c(0,maximum)),
                          angularaxis=list(thetaunit='radians', direction="clockwise", rotation=0)), 
               title=plot_title, legend=list(orientation="h", xanchor="center", x=0.5))
      Sys.sleep(0)
      print(p)
    }else if (exists("sig_diff2")){
      p=plot_ly(type='scatterpolar', mode='lines') %>% 
        add_trace(theta=p1$fv$theta_uncut_z*180/pi, r=p1$fv$fit, line=list(color="blue", width=2.5), name=names_smooths[1]) %>% 
        add_trace(theta=p1$fv$theta_uncut_z*180/pi, r=p1$fv$ul, line=list(color="blue", dash="dot", width=0.5), showlegend=FALSE) %>% 
        add_trace(theta=p1$fv$theta_uncut_z*180/pi, r=p1$fv$ll, line=list(color="blue", dash="dot", width=0.5), showlegend=FALSE) %>% 
        add_trace(theta=p2$fv$theta_uncut_z*180/pi, r=p2$fv$fit, line=list(color="red", width=2.5), name=names_smooths[2]) %>% 
        add_trace(theta=p2$fv$theta_uncut_z*180/pi, r=p2$fv$ul, line=list(color="red", dash="dot", width=0.5), showlegend=FALSE) %>% 
        add_trace(theta=p2$fv$theta_uncut_z*180/pi, r=p2$fv$ll, line=list(color="red", dash="dot", width=0.5), showlegend=FALSE) %>%
        add_trace(theta=seq(sig_diff1[1]*180/pi, sig_diff1[2]*180/pi, length.out=20), 
                  r=c(0, rep(maximum, 18), 0), line=list(color="black", width=0.5), fill="toself", fillcolor=rgb(0,0,0,max=255,alpha=25), showlegend=FALSE) %>%
        add_trace(theta=seq(sig_diff2[1]*180/pi, sig_diff2[2]*180/pi, length.out=20), 
                  r=c(0, rep(maximum, 18), 0), line=list(color="black", width=0.5), fill="toself", fillcolor=rgb(0,0,0,max=255,alpha=25), showlegend=FALSE) %>%
        layout(polar=list(sector=c(20,160), radialaxis=list(angle=90, range=c(0,maximum)),
                          angularaxis=list(thetaunit='radians', direction="clockwise", rotation=0)), 
               title=plot_title, legend=list(orientation="h", xanchor="center", x=0.5))
      Sys.sleep(0)
      print(p)
    }else{
      p=plot_ly(type='scatterpolar', mode='lines') %>% 
        add_trace(theta=p1$fv$theta_uncut_z*180/pi, r=p1$fv$fit, line=list(color="blue", width=2.5), name=names_smooths[1]) %>% 
        add_trace(theta=p1$fv$theta_uncut_z*180/pi, r=p1$fv$ul, line=list(color="blue", dash="dot", width=0.5), showlegend=FALSE) %>% 
        add_trace(theta=p1$fv$theta_uncut_z*180/pi, r=p1$fv$ll, line=list(color="blue", dash="dot", width=0.5), showlegend=FALSE) %>% 
        add_trace(theta=p2$fv$theta_uncut_z*180/pi, r=p2$fv$fit, line=list(color="red", width=2.5), name=names_smooths[2]) %>% 
        add_trace(theta=p2$fv$theta_uncut_z*180/pi, r=p2$fv$ul, line=list(color="red", dash="dot", width=0.5), showlegend=FALSE) %>% 
        add_trace(theta=p2$fv$theta_uncut_z*180/pi, r=p2$fv$ll, line=list(color="red", dash="dot", width=0.5), showlegend=FALSE) %>%
        add_trace(theta=seq(sig_diff1[1]*180/pi, sig_diff1[2]*180/pi, length.out=20), 
                  r=c(0, rep(maximum, 18), 0), line=list(color="black", width=0.5), fill="toself", fillcolor=rgb(0,0,0,max=255,alpha=25), showlegend=FALSE) %>%
        layout(polar=list(sector=c(20,160), radialaxis=list(angle=90, range=c(0,maximum)),
                          angularaxis=list(thetaunit='radians', direction="clockwise", rotation=0)), 
               title=plot_title, legend=list(orientation="h", xanchor="center", x=0.5))
      Sys.sleep(0)
      print(p)
    }
  }
}

2 load in data

df <- read.csv("all_data_NZE_Tongan_w_context_checked_7_March_2019_not_cut_checked_27_Apr.csv", sep=',', stringsAsFactors = F)
# remove empty column
df$X = NULL
df$tokenPooled <- factor(df$tokenPooled)
df$subject <- factor(df$subject)
df$sex <- factor(df$sex)
df$native_lg <- factor(df$native_lg)
# updated 
df$playing_proficiency[df$playing_proficiency == "intermediate"] <- "amateur"
df$playing_proficiency <- factor(df$playing_proficiency, levels = c("amateur","semi-professional","professional"))
df$block <- factor(df$block)
df$theta_uncut_z <- as.numeric(df$theta_uncut_z)
# updated
df$note_intensity <- factor(df$note_intensity, levels = c("piano","mezzopiano","mezzoforte","forte"))
df$tokenPooled <- factor(df$tokenPooled)

3 predict average values for note productions NZE vs Tongan

3.1 Data manipulation

dfSummary <- group_by(df, subject, playing_proficiency, activity, native_lg, tokenPooled, theta_uncut_z_group = cut(theta_uncut_z,breaks=100)) %>% summarise(rhoVar = sd(rho_uncut_z,na.rm=TRUE))
dfSummary$theta_uncut_z = as.character(dfSummary$theta_uncut_z_group)
dfSummary$theta_uncut_z = gsub("\\[|\\]|\\(|\\)", "",dfSummary$theta_uncut_z)
dfSummary$theta_uncut_z = strsplit(dfSummary$theta_uncut_z,",")
dfSummary$theta_uncut_z2 = 1
for(i in c(1:nrow(dfSummary)))
{
  dfSummary$theta_uncut_z2[i] = mean(as.numeric(unlist(dfSummary$theta_uncut_z[i])))
}
dfSummary$theta_uncut_z=dfSummary$theta_uncut_z2
dfSummary$playing_proficiency.ord <- as.ordered(dfSummary$playing_proficiency)
contrasts(dfSummary$playing_proficiency.ord) <- "contr.treatment"
dfSummary$activity.ord <- as.ordered(dfSummary$activity)
contrasts(dfSummary$activity.ord) <- "contr.treatment"
dfSummary$langNoteInt <- interaction(dfSummary$native_lg, 
dfSummary$tokenPooled)
dfSummary$langNoteInt.ord <- as.ordered(dfSummary$langNoteInt)
contrasts(dfSummary$langNoteInt.ord) <- "contr.treatment"
dfSummary$native_lg.ord <- as.ordered(dfSummary$native_lg)
contrasts(dfSummary$native_lg.ord) <- "contr.treatment"
dfSummary$tokenPooled.ord <- as.ordered(dfSummary$tokenPooled)
contrasts(dfSummary$tokenPooled.ord) <- "contr.treatment"
dfSummary = na.omit(dfSummary)
for(i in unique(dfSummary$subject))
{
  dfSubject = subset(dfSummary,dfSummary$subject == i)
  for(j in unique(dfSubject$tokenPooled))
  {
    dfSummary$start[dfSummary$subject == i & dfSummary$tokenPooled == j] <- 
      dfSummary$theta_uncut_z[dfSummary$subject == i & dfSummary$tokenPooled == j] ==
      min(dfSummary$theta_uncut_z[dfSummary$subject == i & dfSummary$tokenPooled == j])
  }
}

3.2 Model specification

if (run_models == TRUE){
  mdl.sys.time1 <- system.time(VAR.gam.noAR.Mod1 <- bam(rhoVar ~ langNoteInt.ord + s(theta_uncut_z, bs="cr", k=10) + s(theta_uncut_z, bs="cr", k=10, by=langNoteInt.ord) +s(theta_uncut_z, subject, bs="fs", k=10, m=1, by=langNoteInt.ord), data=dfSummary, discrete=TRUE, nthreads=ncores, method = "fREML"))
   mdl.sys.time1
   # save model & model summary so they can be reloaded later
   saveRDS(VAR.gam.noAR.Mod1, paste0(output_dir,"/VAR.gam.noAR.Mod1.rds"))
   capture.output(summary(VAR.gam.noAR.Mod1), 
                  file = paste0(output_dir,"/summary_VAR.gam.noAR.Mod1.txt"))
}else{
  # reload model from output_dir
  VAR.gam.noAR.Mod1 = readRDS(paste0(output_dir,"/VAR.gam.noAR.Mod1.rds"))
}
if (run_models == TRUE){
  # set start_value_rho from VAR.gam.noAR.Mod1 as rho_est for the next model
  rho_est <- start_value_rho(VAR.gam.noAR.Mod1)
  mdl.sys.time2 <- system.time(VAR.gam.AR.Mod2 <- bam(rhoVar ~ langNoteInt.ord + s(theta_uncut_z, bs="cr", k=10) + s(theta_uncut_z, bs="cr", k=10, by=langNoteInt.ord) +  s(theta_uncut_z, subject, bs="fs", k=10, m=1, by=langNoteInt.ord), AR.start=dfSummary$start, rho=rho_est, data=dfSummary, discrete=TRUE, nthreads=ncores, method = "fREML"))
  mdl.sys.time2
  # save model & model summary so they can be reloaded later
  saveRDS(VAR.gam.AR.Mod2, paste0(output_dir,"/VAR.gam.AR.Mod2.rds"))
  capture.output(summary(VAR.gam.AR.Mod2), 
                 file = paste0(output_dir,"/summary_VAR.gam.AR.Mod2.txt"))
}else{
  # reload model from output_dir
  VAR.gam.AR.Mod2 = readRDS(paste0(output_dir,"/VAR.gam.AR.Mod2.rds"))
}

3.3 Model inspection

3.3.1 Checking ACF

3.3.1.1 Full

acf_resid(VAR.gam.AR.Mod2, main = "Average ACF AR", cex.lab=1.5, cex.axis=1.5)

3.3.1.2 By Token

acf_resid(VAR.gam.AR.Mod2,split_pred=list(dfSummary$tokenPooled),main = "Average ACF AR by tokenPooled",cex.lab=1.5,cex.axis=1.5)

3.3.1.3 By Language

acf_resid(VAR.gam.AR.Mod2,split_pred=list(dfSummary$native_lg),main = "Average ACF AR by native language",cex.lab=1.5,cex.axis=1.5)

3.3.2 Summary of the model

summary(VAR.gam.AR.Mod2)

Family: gaussian 
Link function: identity 

Formula:
rhoVar ~ langNoteInt.ord + s(theta_uncut_z, bs = "cr", 
    k = 10) + s(theta_uncut_z, bs = "cr", k = 10, by = langNoteInt.ord) + 
    s(theta_uncut_z, subject, bs = "fs", k = 10, m = 1, 
        by = langNoteInt.ord)

Parametric coefficients:
                           Estimate Std. Error t value Pr(>|t|)    
(Intercept)                 8.00231    0.10777  74.256  < 2e-16 ***
langNoteInt.ordNZE.É\u0090      -0.97363    1.61302  -0.604 0.546111    
langNoteInt.ordNZE.É\u0090Ë\u0090    -2.64035    1.85143  -1.426 0.153854    
langNoteInt.ordNZE.É’      -2.46702    2.76454  -0.892 0.372202    
langNoteInt.ordNZE.Bb2     -1.98406    1.29226  -1.535 0.124720    
langNoteInt.ordTongan.Bb2  -1.41836    1.41220  -1.004 0.315220    
langNoteInt.ordNZE.Bb3     -2.98062    1.73953  -1.713 0.086644 .  
langNoteInt.ordTongan.Bb3  -3.58208    2.20014  -1.628 0.103519    
langNoteInt.ordNZE.D4      -4.14783    1.66244  -2.495 0.012604 *  
langNoteInt.ordTongan.D4   -5.16170    1.97641  -2.612 0.009018 ** 
langNoteInt.ordNZE.e        0.50229    1.63619   0.307 0.758855    
langNoteInt.ordTongan.e    -2.15181    0.92909  -2.316 0.020568 *  
langNoteInt.ordNZE.É™      -0.07156    1.26678  -0.056 0.954950    
langNoteInt.ordNZE.É™#     -1.81077    0.40448  -4.477 7.63e-06 ***
langNoteInt.ordNZE.É›      -1.25424    1.83884  -0.682 0.495196    
langNoteInt.ordNZE.ɘ       0.04167    2.05228   0.020 0.983800    
langNoteInt.ordNZE.F3      -3.09511    1.73763  -1.781 0.074894 .  
langNoteInt.ordTongan.F3   -2.53651    1.38509  -1.831 0.067074 .  
langNoteInt.ordNZE.F4      -5.29805    1.44440  -3.668 0.000245 ***
langNoteInt.ordTongan.F4   -8.89295    2.81838  -3.155 0.001606 ** 
langNoteInt.ordTongan.i    -1.95829    1.33837  -1.463 0.143434    
langNoteInt.ordNZE.iË\u0090    -29.28668   16.83868  -1.739 0.082009 .  
langNoteInt.ordTongan.o    -0.92331    5.02574  -0.184 0.854239    
langNoteInt.ordNZE.oË\u0090     -2.13673    1.64411  -1.300 0.193747    
langNoteInt.ordNZE.ɵË\u0090    -1.62651    2.37389  -0.685 0.493247    
langNoteInt.ordTongan.u    -6.83503   10.71965  -0.638 0.523732    
langNoteInt.ordNZE.ʉË\u0090    -3.21586    3.80494  -0.845 0.398022    
langNoteInt.ordNZE.ÊŠ      -1.02949    3.74579  -0.275 0.783444    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
                                                      edf  Ref.df
s(theta_uncut_z)                                    8.872   8.936
s(theta_uncut_z):langNoteInt.ordNZE.É\u0090              5.546   5.946
s(theta_uncut_z):langNoteInt.ordNZE.É\u0090Ë\u0090            5.237   5.619
s(theta_uncut_z):langNoteInt.ordNZE.É’              4.973   5.297
s(theta_uncut_z):langNoteInt.ordNZE.Bb2             4.453   4.756
s(theta_uncut_z):langNoteInt.ordTongan.Bb2          5.678   6.024
s(theta_uncut_z):langNoteInt.ordNZE.Bb3             1.003   1.004
s(theta_uncut_z):langNoteInt.ordTongan.Bb3          7.791   8.078
s(theta_uncut_z):langNoteInt.ordNZE.D4              2.924   3.092
s(theta_uncut_z):langNoteInt.ordTongan.D4           3.844   4.110
s(theta_uncut_z):langNoteInt.ordNZE.e               6.911   7.197
s(theta_uncut_z):langNoteInt.ordTongan.e            6.022   6.443
s(theta_uncut_z):langNoteInt.ordNZE.É™              8.038   8.378
s(theta_uncut_z):langNoteInt.ordNZE.É™#             8.412   8.698
s(theta_uncut_z):langNoteInt.ordNZE.É›              6.415   6.796
s(theta_uncut_z):langNoteInt.ordNZE.ɘ              7.059   7.409
s(theta_uncut_z):langNoteInt.ordNZE.F3              1.001   1.001
s(theta_uncut_z):langNoteInt.ordTongan.F3           6.412   6.792
s(theta_uncut_z):langNoteInt.ordNZE.F4              4.088   4.346
s(theta_uncut_z):langNoteInt.ordTongan.F4           4.886   5.195
s(theta_uncut_z):langNoteInt.ordTongan.i            7.799   8.150
s(theta_uncut_z):langNoteInt.ordNZE.iË\u0090             6.118   6.552
s(theta_uncut_z):langNoteInt.ordTongan.o            6.065   6.458
s(theta_uncut_z):langNoteInt.ordNZE.oË\u0090             5.824   6.194
s(theta_uncut_z):langNoteInt.ordNZE.ɵË\u0090            1.001   1.001
s(theta_uncut_z):langNoteInt.ordTongan.u            6.249   6.720
s(theta_uncut_z):langNoteInt.ordNZE.ʉË\u0090            1.001   1.001
s(theta_uncut_z):langNoteInt.ordNZE.ÊŠ              2.316   2.436
s(theta_uncut_z,subject):langNoteInt.ordNZE.É\u0090     70.857 100.000
s(theta_uncut_z,subject):langNoteInt.ordNZE.É\u0090Ë\u0090   70.884 100.000
s(theta_uncut_z,subject):langNoteInt.ordNZE.É’     79.953 100.000
s(theta_uncut_z,subject):langNoteInt.ordNZE.Bb2    66.321  99.000
s(theta_uncut_z,subject):langNoteInt.ordTongan.Bb2 66.393  99.000
s(theta_uncut_z,subject):langNoteInt.ordNZE.Bb3    71.568  99.000
s(theta_uncut_z,subject):langNoteInt.ordTongan.Bb3 69.313  99.000
s(theta_uncut_z,subject):langNoteInt.ordNZE.D4     67.242  99.000
s(theta_uncut_z,subject):langNoteInt.ordTongan.D4  68.879  99.000
s(theta_uncut_z,subject):langNoteInt.ordNZE.e      66.160  99.000
s(theta_uncut_z,subject):langNoteInt.ordTongan.e   61.414  99.000
s(theta_uncut_z,subject):langNoteInt.ordNZE.É™     59.777  99.000
s(theta_uncut_z,subject):langNoteInt.ordNZE.É™#    48.910  99.000
s(theta_uncut_z,subject):langNoteInt.ordNZE.É›     69.393  99.000
s(theta_uncut_z,subject):langNoteInt.ordNZE.ɘ     71.632  99.000
s(theta_uncut_z,subject):langNoteInt.ordNZE.F3     75.905  99.000
s(theta_uncut_z,subject):langNoteInt.ordTongan.F3  66.936  99.000
s(theta_uncut_z,subject):langNoteInt.ordNZE.F4     61.214  99.000
s(theta_uncut_z,subject):langNoteInt.ordTongan.F4  66.915  99.000
s(theta_uncut_z,subject):langNoteInt.ordTongan.i   62.482  99.000
s(theta_uncut_z,subject):langNoteInt.ordNZE.iË\u0090    79.667  99.000
s(theta_uncut_z,subject):langNoteInt.ordTongan.o   82.606 100.000
s(theta_uncut_z,subject):langNoteInt.ordNZE.oË\u0090    68.673  99.000
s(theta_uncut_z,subject):langNoteInt.ordNZE.ɵË\u0090   77.097  99.000
s(theta_uncut_z,subject):langNoteInt.ordTongan.u   84.433  99.000
s(theta_uncut_z,subject):langNoteInt.ordNZE.ʉË\u0090   81.023  99.000
s(theta_uncut_z,subject):langNoteInt.ordNZE.ÊŠ     77.375  99.000
                                                        F  p-value
s(theta_uncut_z)                                   70.489  < 2e-16
s(theta_uncut_z):langNoteInt.ordNZE.É\u0090              1.445 0.169930
s(theta_uncut_z):langNoteInt.ordNZE.É\u0090Ë\u0090            2.493 0.025088
s(theta_uncut_z):langNoteInt.ordNZE.É’              1.250 0.205920
s(theta_uncut_z):langNoteInt.ordNZE.Bb2             1.584 0.301922
s(theta_uncut_z):langNoteInt.ordTongan.Bb2          2.210 0.041133
s(theta_uncut_z):langNoteInt.ordNZE.Bb3             0.093 0.762294
s(theta_uncut_z):langNoteInt.ordTongan.Bb3          3.113 0.000804
s(theta_uncut_z):langNoteInt.ordNZE.D4              1.113 0.345495
s(theta_uncut_z):langNoteInt.ordTongan.D4           1.625 0.157326
s(theta_uncut_z):langNoteInt.ordNZE.e               3.014 0.004660
s(theta_uncut_z):langNoteInt.ordTongan.e            1.728 0.073853
s(theta_uncut_z):langNoteInt.ordNZE.É™              8.910 1.75e-11
s(theta_uncut_z):langNoteInt.ordNZE.É™#            18.668  < 2e-16
s(theta_uncut_z):langNoteInt.ordNZE.É›              2.142 0.067326
s(theta_uncut_z):langNoteInt.ordNZE.ɘ              3.775 0.000484
s(theta_uncut_z):langNoteInt.ordNZE.F3              0.057 0.812421
s(theta_uncut_z):langNoteInt.ordTongan.F3           2.265 0.018740
s(theta_uncut_z):langNoteInt.ordNZE.F4              1.708 0.170357
s(theta_uncut_z):langNoteInt.ordTongan.F4           2.669 0.033058
s(theta_uncut_z):langNoteInt.ordTongan.i            3.002 0.001150
s(theta_uncut_z):langNoteInt.ordNZE.iË\u0090             2.756 0.006887
s(theta_uncut_z):langNoteInt.ordTongan.o            2.203 0.048544
s(theta_uncut_z):langNoteInt.ordNZE.oË\u0090             2.666 0.012507
s(theta_uncut_z):langNoteInt.ordNZE.ɵË\u0090            0.026 0.873668
s(theta_uncut_z):langNoteInt.ordTongan.u            2.042 0.055998
s(theta_uncut_z):langNoteInt.ordNZE.ʉË\u0090            0.413 0.520613
s(theta_uncut_z):langNoteInt.ordNZE.ÊŠ              0.180 0.783067
s(theta_uncut_z,subject):langNoteInt.ordNZE.É\u0090      6.865  < 2e-16
s(theta_uncut_z,subject):langNoteInt.ordNZE.É\u0090Ë\u0090    7.818  < 2e-16
s(theta_uncut_z,subject):langNoteInt.ordNZE.É’     18.664  < 2e-16
s(theta_uncut_z,subject):langNoteInt.ordNZE.Bb2    10.303  < 2e-16
s(theta_uncut_z,subject):langNoteInt.ordTongan.Bb2 18.281  < 2e-16
s(theta_uncut_z,subject):langNoteInt.ordNZE.Bb3    14.466  < 2e-16
s(theta_uncut_z,subject):langNoteInt.ordTongan.Bb3 22.632  < 2e-16
s(theta_uncut_z,subject):langNoteInt.ordNZE.D4     13.933  < 2e-16
s(theta_uncut_z,subject):langNoteInt.ordTongan.D4  25.480  < 2e-16
s(theta_uncut_z,subject):langNoteInt.ordNZE.e      12.291  < 2e-16
s(theta_uncut_z,subject):langNoteInt.ordTongan.e    8.573  < 2e-16
s(theta_uncut_z,subject):langNoteInt.ordNZE.É™     14.365  < 2e-16
s(theta_uncut_z,subject):langNoteInt.ordNZE.É™#     3.117  < 2e-16
s(theta_uncut_z,subject):langNoteInt.ordNZE.É›      8.070  < 2e-16
s(theta_uncut_z,subject):langNoteInt.ordNZE.ɘ     13.504  < 2e-16
s(theta_uncut_z,subject):langNoteInt.ordNZE.F3     15.014  < 2e-16
s(theta_uncut_z,subject):langNoteInt.ordTongan.F3  20.935  < 2e-16
s(theta_uncut_z,subject):langNoteInt.ordNZE.F4     11.949  < 2e-16
s(theta_uncut_z,subject):langNoteInt.ordTongan.F4  26.468  < 2e-16
s(theta_uncut_z,subject):langNoteInt.ordTongan.i   14.004  < 2e-16
s(theta_uncut_z,subject):langNoteInt.ordNZE.iË\u0090    45.669  < 2e-16
s(theta_uncut_z,subject):langNoteInt.ordTongan.o   22.896  < 2e-16
s(theta_uncut_z,subject):langNoteInt.ordNZE.oË\u0090    13.644  < 2e-16
s(theta_uncut_z,subject):langNoteInt.ordNZE.ɵË\u0090   12.675  < 2e-16
s(theta_uncut_z,subject):langNoteInt.ordTongan.u   43.466  < 2e-16
s(theta_uncut_z,subject):langNoteInt.ordNZE.ʉË\u0090   22.778  < 2e-16
s(theta_uncut_z,subject):langNoteInt.ordNZE.ÊŠ     28.619  < 2e-16
                                                      
s(theta_uncut_z)                                   ***
s(theta_uncut_z):langNoteInt.ordNZE.É\u0090                
s(theta_uncut_z):langNoteInt.ordNZE.É\u0090Ë\u0090           *  
s(theta_uncut_z):langNoteInt.ordNZE.É’                
s(theta_uncut_z):langNoteInt.ordNZE.Bb2               
s(theta_uncut_z):langNoteInt.ordTongan.Bb2         *  
s(theta_uncut_z):langNoteInt.ordNZE.Bb3               
s(theta_uncut_z):langNoteInt.ordTongan.Bb3         ***
s(theta_uncut_z):langNoteInt.ordNZE.D4                
s(theta_uncut_z):langNoteInt.ordTongan.D4             
s(theta_uncut_z):langNoteInt.ordNZE.e              ** 
s(theta_uncut_z):langNoteInt.ordTongan.e           .  
s(theta_uncut_z):langNoteInt.ordNZE.É™             ***
s(theta_uncut_z):langNoteInt.ordNZE.É™#            ***
s(theta_uncut_z):langNoteInt.ordNZE.É›             .  
s(theta_uncut_z):langNoteInt.ordNZE.ɘ             ***
s(theta_uncut_z):langNoteInt.ordNZE.F3                
s(theta_uncut_z):langNoteInt.ordTongan.F3          *  
s(theta_uncut_z):langNoteInt.ordNZE.F4                
s(theta_uncut_z):langNoteInt.ordTongan.F4          *  
s(theta_uncut_z):langNoteInt.ordTongan.i           ** 
s(theta_uncut_z):langNoteInt.ordNZE.iË\u0090            ** 
s(theta_uncut_z):langNoteInt.ordTongan.o           *  
s(theta_uncut_z):langNoteInt.ordNZE.oË\u0090            *  
s(theta_uncut_z):langNoteInt.ordNZE.ɵË\u0090              
s(theta_uncut_z):langNoteInt.ordTongan.u           .  
s(theta_uncut_z):langNoteInt.ordNZE.ʉË\u0090              
s(theta_uncut_z):langNoteInt.ordNZE.ÊŠ                
s(theta_uncut_z,subject):langNoteInt.ordNZE.É\u0090     ***
s(theta_uncut_z,subject):langNoteInt.ordNZE.É\u0090Ë\u0090   ***
s(theta_uncut_z,subject):langNoteInt.ordNZE.É’     ***
s(theta_uncut_z,subject):langNoteInt.ordNZE.Bb2    ***
s(theta_uncut_z,subject):langNoteInt.ordTongan.Bb2 ***
s(theta_uncut_z,subject):langNoteInt.ordNZE.Bb3    ***
s(theta_uncut_z,subject):langNoteInt.ordTongan.Bb3 ***
s(theta_uncut_z,subject):langNoteInt.ordNZE.D4     ***
s(theta_uncut_z,subject):langNoteInt.ordTongan.D4  ***
s(theta_uncut_z,subject):langNoteInt.ordNZE.e      ***
s(theta_uncut_z,subject):langNoteInt.ordTongan.e   ***
s(theta_uncut_z,subject):langNoteInt.ordNZE.É™     ***
s(theta_uncut_z,subject):langNoteInt.ordNZE.É™#    ***
s(theta_uncut_z,subject):langNoteInt.ordNZE.É›     ***
s(theta_uncut_z,subject):langNoteInt.ordNZE.ɘ     ***
s(theta_uncut_z,subject):langNoteInt.ordNZE.F3     ***
s(theta_uncut_z,subject):langNoteInt.ordTongan.F3  ***
s(theta_uncut_z,subject):langNoteInt.ordNZE.F4     ***
s(theta_uncut_z,subject):langNoteInt.ordTongan.F4  ***
s(theta_uncut_z,subject):langNoteInt.ordTongan.i   ***
s(theta_uncut_z,subject):langNoteInt.ordNZE.iË\u0090    ***
s(theta_uncut_z,subject):langNoteInt.ordTongan.o   ***
s(theta_uncut_z,subject):langNoteInt.ordNZE.oË\u0090    ***
s(theta_uncut_z,subject):langNoteInt.ordNZE.ɵË\u0090   ***
s(theta_uncut_z,subject):langNoteInt.ordTongan.u   ***
s(theta_uncut_z,subject):langNoteInt.ordNZE.ʉË\u0090   ***
s(theta_uncut_z,subject):langNoteInt.ordNZE.ÊŠ     ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.868   Deviance explained = 88.2%
fREML =  36912  Scale est. = 2.0136    n = 19150

3.3.3 Plotting model predictions

3.3.3.1 Overall variance differences

And finally overall variance differences by language

# make copy of dfNotes
dat1 = dfSummary
dat1$predicted_values = predict(VAR.gam.AR.Mod2)
# plot in polar coordinates using plotly
dat1_NZE = dat1[dat1$native_lg == "NZE",]
dat1_Tongan = dat1[dat1$native_lg == "Tongan",]
# estimate smooths using R's generic predict.smooth.spline function
smooth_NZE=list(theta=seq(min(dat1_NZE$theta_uncut_z)*180/pi, max(dat1_NZE$theta_uncut_z)*180/pi, length=100), r=predict(smooth.spline(dat1_NZE$theta_uncut_z, dat1_NZE$predicted_values),
 seq(min(dat1_NZE$theta_uncut_z), max(dat1_NZE$theta_uncut_z), length=100))$y, line=list(color="blue", dash="dash"))
smooth_Tongan=list(theta=seq(min(dat1_Tongan$theta_uncut_z)*180/pi, max(dat1_Tongan$theta_uncut_z)*180/pi, length=100), r=predict(smooth.spline(dat1_Tongan$theta_uncut_z, dat1_Tongan$predicted_values),
 seq(min(dat1_Tongan$theta_uncut_z), max(dat1_Tongan$theta_uncut_z), length=100))$y, line=list(color="red", dash=""))
# set Rho max to the max of the predicted_values + 5
max_predictions = max(max(smooth_NZE$r), max(smooth_Tongan$r))
maximum=max_predictions+5
rm(max_predictions)
p = plot_ly(type='scatterpolar', mode='lines') %>%
add_trace(theta=smooth_NZE$theta, r=smooth_NZE$r, line=list(color=smooth_NZE$line$color[[1]], width=2.5, dash=smooth_NZE$line$dash[[1]]), name="overall average of NZE notes") %>%
  add_trace(theta=smooth_Tongan$theta, r=smooth_Tongan$r, line=list(color=smooth_Tongan$line$color[[1]], width=2.5, dash=smooth_Tongan$line$dash[[1]]), name="overall average of Tongan notes") %>%
  layout(polar=list(sector=c(20,160), radialaxis=list(angle=90, range=c(0,maximum)), angularaxis=list(thetaunit='radians', direction="clockwise", rotation=0)), title="Overall average smooths for NZE vs Tongan note productions", legend=list(orientation="h", xanchor="center", x=0.5))
p

rm(dat1, dat1_NZE, dat1_Tongan, maximum, smooth_NZE, smooth_Tongan)